package algorithm.sphericharmonics;

import java.util.Hashtable;

import algorithm.fft.Complex;
import algorithm.fft.FFT;

import edu.emory.mathcs.jtransforms.fft.FloatFFT_1D;

public final class VoxelizedSphere {
	public static final int lmax=17;
	private int sphereIndex; 
	private float radius; 
	private float delta_theta; 
	private float fromAngle=Float.MAX_VALUE; 
	private float toAngle=Float.MIN_VALUE;
	private static final float PI=(float)Math.PI;
	private static final float TWOPI=(float)Math.PI*2;
	public static final int[] lmIndex_offset=new int[lmax]; 
	static {
		int k=0; 
		for (int l=0; l<lmax;l++){
			lmIndex_offset[l]=k; 
			k=k+l+1; 
		}
	}
	private float[] l2Normal=new float[lmax]; 
	
	private Complex[] alm; 
	private int numberOfGrid; 
	private float delta_phi; 
	private final Hashtable<Integer,LongitudeCircle> rows=new Hashtable<Integer,LongitudeCircle>();
	public VoxelizedSphere(float r,int numOfGrid,int sphIndex){
		radius=r;
		this.numberOfGrid=numOfGrid; 
		delta_theta=PI/numberOfGrid;
		delta_phi=delta_theta; 
		this.sphereIndex=sphIndex; 
		alm=new Complex[((lmax+1)*lmax)/2]; 
		for (int i=0; i<alm.length; i++){
			alm[i]=new Complex(0,0); 
		}
	}
	public int numberOfGrid(){
		return numberOfGrid;  
	}
	public int getLngNumOfGrid(){
		return 2*numberOfGrid; 
	}
	public final void addPixel(float theta, float phi){
		if (phi>toAngle){
			toAngle=phi; 
		}
		int index_phi=(int)(phi/this.delta_phi); 
		if (phi<fromAngle){
			fromAngle=phi; 
		}
		LongitudeCircle row=getRow(theta); 
		row.add(index_phi); 
	}
	private Complex alm_negative=new Complex(0,0); 
	private final Complex alm(int l, int m){
		int mm=m>=0?m:-m;
		int index=lmIndex_offset[l]+mm; 
		Complex c=this.alm[index];
		if (m>=0){
			return c;  
		}else {
			if (mm%2==0){// even
				alm_negative.re=c.re; 
				alm_negative.im=-c.im; 
			}else {
				alm_negative.re=-c.re; 
				alm_negative.im=c.im; 				
			}
			return alm_negative; 
		}
	}
	public static final int almIndex(int l,int m) {
		return  lmIndex_offset[l]+m; 
	}
	public Complex a(int l, int m){
		//float range=toAngle-fromAngle+delta_theta/2; 
		//int zero_index=(int)(fromAngle/delta_theta);  
		int zero_index=0; 
		int n=this.getLngNumOfGrid();   
		if (n<=0){
			return new Complex(0,0); // problematic, how come; 
		}  
		int buffer_size=Math.max(n, lmax); 
		float[] buffer=new float[2*buffer_size];
		Complex c=new Complex(0,0);		 
		for (LongitudeCircle circle: rows.values()){
			java.util.Arrays.fill(buffer,0, buffer.length, 0); 
			IntList k=circle.getLngAngles();
			int pos=0; 
			//System.out.println("k.size():"+k.size()+" zero_index:"+zero_index); 
			for (int i=0; i<k.size(); i++){
				pos=k.get(i)-zero_index; 
				buffer[pos]=1;  // basically a boolean array; 
			}
			float theta=circle.getLatAngle(); 
			Complex d=new Complex(0,0); 
			for (int i=0; i<n; i++){ 
				if (buffer[i]>0.5f){
					float phi=i*circle.getDeltaPhi();
					Complex sh=SphericalHarmonics.yfunct(l, m, theta, phi).toComplex();
					d.re=d.re+sh.re; 
					d.im=d.im+sh.im; 
				}
			}
			c.re=c.re+d.re; 
			c.im=c.im+d.im;
			//System.out.println("========"+l+","+m+" circle:"+circle.getIndex()+"============");
			//System.out.println("alm:"+d); 
			//System.out.println("------------------------------------"); 
		}
		return c; 
	}
	public final float[] transform(){ 
		fromAngle=0; 
		//toAngle=(float)Math.PI*2; 		
		//int zero_index=(int)(fromAngle/delta_theta);
		//int end_index=(int)(toAngle/delta_theta); 
		int n=this.getLngNumOfGrid(); 
		
		if (n<lmax){
			throw new RuntimeException("longitute grid "+n+" <lmax "+lmax); 
		} 
		float[] buffer=new float[2*n];
		float[] qValues; 
		float re_start;  
		float im_start;  
		FloatFFT_1D fft_processor=new FloatFFT_1D(n); //here is where the problem is. 
		float[] norm = new float[lmax]; 
		/* calculate sqrt((2*l + 1)/(4*PI)) */
		float temp = (float)Math.sqrt(1/(4*Util.PI));
		for(int l = 0; l < lmax; l++){
		    norm[l] = temp*(float)Math.sqrt(2*l+1);
		}
		float re_m=0;  
		float im_m=0;
		float re; 
		float im; 
		float q_value;
		
		for (LongitudeCircle circle: rows.values()){
			IntList k=circle.getLngAngles();
			int pos=0;
			//System.out.println("theta:"+circle.getLatAngle());
			java.util.Arrays.fill(buffer,0, buffer.length, 0);
			for (int i=0; i<k.size(); i++){
				pos=k.get(i);//-zero_index; 
				buffer[pos]=1;  // basically a boolean array;
			} 
			fft_processor.realForwardFull(buffer); // do the fft to the phi data;
			//fft_processor.complexInverse(buffer, false); 
			for (int m=0; m<lmax; m++){
				//re_start=(float)Math.cos(fromAngle*m); 
				//im_start=(float)Math.sin(fromAngle*m);
				//qValues=circle.getQValues(m);
				qValues=circle.qFunctions[m].qValues;
				float ffted_re=buffer[2*m]; 
				float ffted_im=buffer[2*m+1]; 
				for (int l=m; l<lmax; l++){				
					int index=lmIndex_offset[l]+m; //almIndex(l,m); only for performance reason; 			
					Complex a_postive=alm[index];					
					q_value=qValues[l-m]; //circle.getValue(l, m);
					/*
					re_m=buffer[2*m]*q_value; 
					im_m=buffer[2*m+1]*q_value;
					re=re_m*re_start-im_m*im_start; 
					im=re_m*im_start+re_start*im_m;
					*/
					re=ffted_re*q_value; 
					im=ffted_im*q_value;
					a_postive.re+=re*norm[l]; 
					a_postive.im+=im*norm[l];
				}
			} 
		}
		float l2norm;
		int index;
		int start; 
		Complex c; 
		for (int l=0; l<lmax; l++){				
			start=lmIndex_offset[l];
			c=alm[start];
			l2norm=(float)(c.re*c.re+c.im*c.im);
			l2Normal[l]=l2norm; 
			for (int m=1; m<=l; m++){
				index=start+m; 
				c=this.alm[index];
				l2norm=(float)(c.re*c.re+c.im*c.im);
				l2Normal[l]+=(l2norm+l2norm); 
			}
			
			/*
			for (int m=-l; m<=l; m++){
				c=this.alm(l, m);  
				l2Normal[l]+=(float)(c.re*c.re+c.im*c.im);
				
				Complex ctest=this.a(l, m); 
				//ctest=ctest.minus(c);
				double error=Math.abs(ctest.re-c.re);//+Math.abs(ctest.im-c.im);
				if (error>0.001){
					System.out.println(l+","+m+" error:"+error+" c:"+c+"/ ctest:"+ctest); 
				}
				
			}*/
			
			l2Normal[l]=(float)Math.sqrt(l2Normal[l]); 
			/*
			if (l2Normal[l]>0.00001){
				System.out.println("f("+l+")="+l2Normal[l]);
			}
			*/
						
		}
		return l2Normal; 
	}
	
	public int size(){
		int size=0; 
		for (LongitudeCircle row:rows.values()){
			size+=row.lngAngles.size();
		}
		if (size==0){
			System.out.println("size:0"); 
		}
		return size; 
	}
	private final LongitudeCircle getRow(float theta){
		int theta_index=(int)(theta/delta_theta);
		LongitudeCircle row=rows.get(theta_index); 
		if (row==null){
			row=new LongitudeCircle(theta,2*this.numberOfGrid, theta_index); 
			rows.put(theta_index, row); 
		}
		return row; 
	}
	public float getRadius() {
		return radius;
	}
	public void setRadius(float radius) {
		this.radius = radius;
	}
	public float getFromAngle() {
		return fromAngle;
	}
	public void setFromAngle(float fromAngle) {
		this.fromAngle = fromAngle;
	}
	public int getSphereIndex() {
		return sphereIndex;
	}
	public static void testRotation(){
		float theta=10*(float)Math.PI/180.f;
		float phi=0; 
		VoxelizedSphere ll=new VoxelizedSphere(1.0f,180,0); 
		//LongitudeCircle ll=new LongitudeCircle(theta,180,0);
		for (int l=0; l<=4;l++){
			System.out.println(l+"========");
			float[] vals=new float[180];
			for (int a=0; a<180; a++){
				theta=(float)(Math.PI/180)*a; 
				float v=0; 
				phi=10*(float)Math.PI/180.f;
				for (int m=-l; m<=l; m++){
					Complex sh=SphericalHarmonics.yfunct(l, m, theta, phi).toComplex();
					v+=sh.re*sh.re+sh.im*sh.im; 
				}
				v=(float)Math.sqrt(v); 
				//System.out.println(v); 
				vals[a]=v;  
			}
			float avg=0; 
			for (int i=0; i<180; i++){
				avg+=vals[i]; 
			}
			avg=avg/vals.length;
			float dev=0; 
			for (int i=0; i<180; i++){
				dev+=(vals[i]-avg)*(vals[i]-avg); 
			}
			dev=(float)Math.sqrt(dev/vals.length);
			System.out.println(" avg :"+avg +" dev:"+dev); 
		}
	}
	public static void main(String[] args){
		
		float r=1.0f; 
		float dtheta=0.74331594f; 
		float theta=2.3201294f;
		float phi=2.2299478f;

		int sphereIndex=1;
		testRotation(); 
		/*
		for (int i=0; i<lmax;i++){
			System.out.println("a("+i+",0) start at:"+almIndex_offset[i]); 
		}
		*/
		/*
		VoxelizedSphere v=new VoxelizedSphere(r, dtheta,sphereIndex);
		v.addPixel(theta, phi); 
		theta=(float)Math.PI/10; 
		//phi=dtheta*20; 
		v.addPixel(theta,phi); 
		v.transform();
		*/
	}
}
